Dust cloud formation in stellar environments 



II. Two-dimensional models for structure formation around AGB stars 

P. Woitke 1 ' 2 and G. Niccolini 3 

1 Sterrewacht Leiden, P.O Box 9513, 2300 RA Leiden, The Netherlands 

2 Zentrum fur Astronomie und Astrophysik, TU Berlin, Hardenbergstrafie 36, D- 10623 Berlin 

3 Observatoire de la Cote dAzur, Departement Fresnel UMR 6528, BP 4229, F-06034 Nice Cedex 4, France 



o 
o 

(N 
u 

< 
oo 



> 
O 

in 

o 

o 

Or 
i 1 

o ■ 
a: 



Received 9 February 2004; accepted date 

Abstract. This paper reports on computational evidence for the formation of cloud-like dust structures around 
C-rich AGB stars. This spatio-temporal structure formation process is caused by a radiative/thermal instability of 
dust forming gases as identified by Woitke et al. (2000). Our 2D (axisymmetric) models combine a time-dependent 
description of the dust formation process according to Gail & Sedlmayr (1988) with detailed, frequency-dependent 
continuum radiative transfer by means of a Monte Carlo method (Niccolini et al. 2003) in an otherwise static 
medium (u=0). These models show that the formation of dust behind already condensed regions, which shield 
the stellar radiation field, is strongly favoured. In the shadow of these clouds, the temperature decreases by several 
hundred Kelvin which triggers the subsequent formation of dust and ensures its thermal stability. Considering an 
initially dust-free gas with small density inhomogeneities, we find that finger-like dust structures develop which 
are cooler than the surroundings and point towards the centre of the radiant emission, similar to the "cometary 
knots" observed in planetary nebulae and star formation regions. Compared to a spherical symmetric reference 
model, the clumpy dust distribution has little effect on the spectral energy distribution, but dominates the optical 
appearance in near IR monochromatic images. 

Key words. Instabilities - Radiative transfer - Dust, extinction - Stars: AGB and post- AGB - Circumstellar 
matter - Stars: mass loss 



1. Introduction 

Dusty gases in space are often remarkably inhomogeneous. 
Numerous observations of various dust forming objects 
like the circumstellar environments of AGB and post- AGB 
stars, R Coronae Borealis stars, planetary nebulae, the 
ejecta of novae and supernovae, and even the hot winds 
generated by Wolf-Rayet stars, have shown that the dust 
forming medium usually possesses a clumpy internal struc- 
ture. Reviews of such observations are given by Lopez 
(1999) and Woitke (2001). 

The best-studied object in that respect is probably the 
infrared carbon star IRC+10216. Several infrared speckle 
observations of its innermost dust formation and wind 
acceleration zone show direct evidence for an irregular, 
possibly cloudy dust distribution around this late-type 
AGB star (Weigelt et al. 1998, Haniff & Buscher 1998, 
Tuthill et al. 2000). Further evidence can be deduced from 
long-term JHKL lightcurves. Concerning the carbon star 
II Lup, Feast et al. (2003) argue for a restricted epoch 
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of dust formation in a limited region along the line of 
sight, in order to explain a large-amplitude long-term de- 
crease in J whereas no corresponding increase in K and 
L was detected. Direct IR imaging and multi-wavelength 
IR lightcurves of the oxygen-rich red giant L2 Pup point 
to similar wind asymmetries (Jura et al. 2002). 

Additional hints to inhomogeneities in the environ- 
ments of AGB stars are given by the patchy SiO maser 
spots observed in oxygen-rich Mira variables (e. g. TX 
Cam, Diamond & Kemball 2003). On a larger scale, CO 
rotational emission lines provide evidence for local den- 
sity enhancements in the winds of late-type stars, e. g. 
concerning the carbon star TX Psc (Heske et al. 1989) or 
in high-resolution observations of the detached shell of TT 
Cyg (Olofsson et al. 2000). Observations with higher spa- 
tial resolution, using instruments like the Vlti, Ngst or 
Alma, can be expected to reveal even more details in the 
near future. 

Spherically symmetric dynamical models for AGB 
winds which include a time-dependent treatment of 
dust formation strongly suggest the formation of ra- 
dial dust shells in more or less regular time inter- 



vals (e.g. Winters et al. 2000, bandin & Homer 2003), 
even if the pulsation of the star is neglected (e.g. 
Fleischer et al. 1995). The question arises whether these 
dust shells remain spherical symmetric (as suggested by 
the ID models), or whether they might break apart into 
clouds due to instabilities. We focus in this paper on the 
second possibility and study the two-dimensional time- 
dependent behaviour of the dust/gas mixture just during 
the formation of a new dust shell. 

dumpiness may also play a vital role for the photo- 
chemistry during the AGB — > PPN — > PN transition phase. 
Opaque clumps can suppress the photodissociation of cer- 
tain molecules like Benzene, such that their existence has 
been proposed to be an indicator for inhomogeneities with 
large density contrasts (Redman et al. 2003). 

Regarding other classes of objects, numerous opaque 
structures have been discovered in proto-planetary neb- 
ulae (e. g. the Helix nebula NGC 7293 and the Eskimo 
nebula NGC 2392) as well as in star formation regions 
(e.g. the Orion nebula). These neutral, dense, probably 
dusty regions, designated as "cometary knots", "globules" 
or "proplyds" , are located at the head of radially aligned 
linear structures in these nebulae and are particularly well 
visible in high resolution images of emission line ratios like 
[Oin]/Ha (O'Dell 2000). 

The physical cause of the observed dumpiness is still 
puzzling. In Paper I of this series (Woitke et al. 2000) , we 
have formulated the hypothesis that a radiative/thermal 
instability in dust forming gases can provoke a self- 
organisation of the matter, which is possibly involved in 
the formation of the observed structures. This instabil- 
ity is characterised by a physical control loop between 
the radiative transfer, which determines the temperature 
structure of the medium, and the dust formation, which 
determines its opacity (see Fig. 1 in Paper I). 

This paper reports on computational evidence for this 
hypothesis. In Sect. 2, we outline the concept of our static 
model for the environment of a C-rich AGB star, which 
combines a time-dependent treatment of dust formation 
with two-dimensional radiative transfer. Section 3 shows 
and discusses the results, including the calculated opti- 
cal appearance of a clumpy dust distribution in the spec- 
tral energy distribution and in monochromatic images. In 
Sect. 4, our conclusions are drawn. 

2. The model 

The spatio-temporal evolution of the dust component 
in the circumstellar environment of a C-rich AGB star 
is simulated by means of a 2D (axisymmetric) time- 
dependent model. Our aim is to investigate whether ra- 
diative/thermal instabilities can excite a spatio-temporal 
self-organisation of the dust forming medium. We do not 
intend to model any specific object. Instead, we are inter- 
ested in what kind of spatial structures possibly develop 
and we want to identify which physical control mecha- 
nisms are responsible for the structure formation. 



According to this purpose, the radiative, chemical and 
physical processes directly involved in the instability must 
be modelled as complete and precise as necessary. In our 
two-dimensional approach, a precise radiative transfer re- 
quires large computational efforts (see Sects. I2~B1 and l2~f)l 
for details) and with our current approach, we are already 
at the computational limit of todays super-computer fa- 
cilities. Therefore, the remaining parts of the model must 
be kept comparably simple. For these reasons, we assume 
that the gas/dust mixture is static (v = 0) and focus only 
on the time-dependent dust-chemistry and the radiative 
transfer. Another reason for this concept is that we want 
to study the pure effect of the instability, which might get 
entangled with other dynamical instabilities, if hydrody- 
namics was included. The question whether or not this 
model is realistic depends on the relation between the 
time-scales of the included processes (radiative transfer, 
chemistry and dust formation) and the hydrodynamical 
time-scale (see discussion in Paper I). If the formation of 
a new dust shell occurs on a much smaller time-scale than 
the hydrodynamical time-scale, the assumption of a quasi- 
static medium seems appropriate 1 . Regarding the param- 
eters of our standard model (see Fig. [2J , in particular the 
stellar luminosity i^ = 3000LQ, it is not sure that enough 
radiation pressure can be provided to drive a dust-driven 
outflow, i. e. the hydrodynamical time-scale could be in 
fact quite large. 

Our approach, at least, allows for a first explo- 
ration of the qualitative two-dimensional behaviour of 
the dust component in AGB star winds which, accord- 
ing to our knowledge, is the first multi-dimensional and 
time-dependent (though static) model for dust formation 
around a cool star which includes the important feedback 
of dust on the temperature structure of the medium via 
radiative transfer effects 2 . 

2.1. Gas density distribution 

The time-dependent simulations of dust formation and 
radiative transfer are carried out for a restricted model 
volume, which covers the outer atmospheric layers of the 
star mainly responsible for the dust formation. Within this 
volume, the hydrogen nuclei particle density [cm -3 ] 
(mass density p = 1.4amu-n{ H j for solar element abun- 
dances) is prescribed and left unchanged during the time- 
dependent calculations, according to our assumption of a 
static medium. We consider an exponentially decreasing 



In the periodically shocked extended atmospheres of long- 
period variables, the hydrodynamical time-scale is of the order 
of one pulsation period. 

2 Other authors have mainly focused on the dynamics 
and the wind acceleration by means of spherically sym- 
metric models (e. g. Winters et al. 2000, Homer et al. 2003, 
Sandin & Homer 2003). Freytag & Homer (2003) have pas- 
sively integrated the dust moment equations in first 3D simu- 
lations, without feedbacks. 
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n (n)(r) = "<h}M -exp ( — h<5(r) 



(1) 



where H p is the density scale height and n(n) (fo) the mean 
hydrogen nuclei particle density at the inner radial bound- 
ary of the model. 5(r) denotes a small spatial density per- 
turbation of the order of 3% (see App. A for details). 

2.2. Chemistry and dust formation 

The dust formation process is simulated as function of 
time t and space r. The dust component is described by 
moments of the size distribution function f(a,r,t), intro- 
duced by Gail & Sedlmayr (1988), which are defined by 



K 3 (r,t) 



OO 

J f(a,r,t) da . 



(2) 



The dust particles are assumed to be spheres of radius a, 
to be composed of a unique solid material (amorphous car- 
bon) and to have a size- independent temperature ?d(r, i). 
ao is the hypothetical monomer radius for graphite and 
at a lower integration boundary in size space, chosen to 
be 10 • ao . The "~ marks quantities defined per H-atom, for 
example / = //n(H> [cm -1 ]. 

The temporal evolution of the dimensionless dust 
moments Kj is described by a system of differential 
equations in conservation form (Gail & Sedlmayr 1988), 
which accounts for nucleation and growth. The method 
has been extended to include thermal evaporation by 
Gauger et al. (1990) 
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J(at) [cm _1 s _1 ] is the flux of dust particles in size space 
through the lower integration boundary at, and r gr the 
growth time scale defined by Eq. (0- For net growth 
(r" 1 > 0), J(at) can be identified with the nucleation 
rate J* (Eq.JHJl. For net evaporation (t -1 < 0), J (at) is 
calculated from the actual size distribution function at at 
and the shrink rate da/dt fEq.H4|). 



3 Our motivation for this type of radial density-dependence 
originates from ID models (e.g. Fleischer et al. 1995), which 
show an exponential rather than a l/r 2 -dependence close to 
the star, just before the formation of a new dust shell. Note, 
that the medium is not in hydrostatic equilibrium, because H p 
is considered as a free parameter. 
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The dust moment equations (Eq.[3J) are completed by 
the carbon element conservation equation. Since we ne- 
glect relative velocities between the gas and dust particles 
(drift) in our model, it is possible to express the element 
conservation by a simple algebraic equation (Eq.EJ, where 
ec is the carbon abundance relative to hydrogen in the gas 
phase, and its (initial) value in the uncondensed case. 
The degree of condensation is defined by 



/cond Q 



(6) 



where eo is the oxygen element abundance. For the calcu- 
lation of the various chemical and surface reaction rates 
involved in the nucleation, growth and evaporation of the 
dust particles, several particle densities in the gas phase 
rife [cm -3 ] are required (see Table ^) , in general all par- 
ticle densities for atoms, electrons, ions and molecules. 
These particle densities are calculated by assuming chem- 
ical equilibrium in the gas phase, which can formally be 
written as 

n k = n k (n {H) ,T g ,e c ) , (7) 

where T g is the gas temperature. Elements other than 
carbon are assumed to have solar abundances (Anders 
& Grevesse 1989). We use 13 elements (H, He, C, N, 
O, Si, Mg, Al, Fe, S, Na, K, Ti) and 150 molecules in 
our chemical equilibrium code with equilibrium constants 
newly fitted to the electronic data of the Janaf tables 
(Chase et al. 1985). The total carbon abundance Cq is a 
free parameter. For the model discussed in this paper, the 
carbon-to-oxygen ratio is assumed to be 2, i. e. = 2 • eo- 
For the calculation of the nucleation rate J*, we as- 
sume homogeneous nucleation of pure carbon clusters C n 
and apply modified classical nucleation theory accord- 
ing to the scheme of Gail et al. (1984). All molecules and 
clusters are assumed to be in LTE with the gas phase. 
Consequently, the nucleation rate is calculated according 
to T g . Besides the directly involved small carbon clusters 
C, C2 and C3, the nucleation rate depends on the concen- 
tration of some additional small hydro-carbon-molecules 



4 A more reliable description of the nucleation process un- 
der astrophysical conditions is still not at hand (see e. g. dis- 
cussion in Andersen et al. 2003). Detailed chemical rate net- 
works, e.g. for PAH formation (Cherchneff et al. 1992), tend 
to predict much to low formation yields. Future progress in 
this field can be achieved by determining individual cluster 
data by ab-initio quantum mechanical calculations, e. g. for 
(Ti0 2 ) x (Jeong 2000) and Al x O y (Patzer et al. 2003). 



which contribute to the growth ot the carbon clusters 

J* = J*(r g ,n c ,7ic 2 ,nc 2 H,wc 2 H 2 ,wc 3 ,«C3H) • (8) 

The net growth of all dust particles included in the dust 
moments (a > a£) is described by a size- independent 
growth time-scale 



./ 



gr (^gj ^d) n c,nc 2 j n C 2 H, ™C 2 H 2 , n C3 , «c 3 h) 
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n r t£ e m r a r I 1 — 



(9) 



n r is the density of the impinging gas particle which initi- 
ates the reaction r, a r is the sticking coefficient and m r the 
number of monomers (carbon atoms) added to the solid 
surface per reaction. For further definitions and expla- 
nations please consult Gail et al. (1984), Gail & Sedlmayr 
(1988) and Gauger et al. (1990). 

The latter factor in Eq. (JjJJ accounts for the reverse 
(evaporation) reactions by applying Milne relations with 
respect to detailed balance 



,rr \ > whcrc 

p vap (T) = 10 6 dyn/cm 2 .exp(^|^ , 

AG(T) = +1.01428- 10 6 /T - 7.23043 • 10 5 

+ 1.63039 • 10 2 • T - 1.75890 • 10~ 3 • T 2 
+ 9.97416- 10~ 8 • T 3 . 



(10) 
(11) 

(12) 



S is the supersaturation ratio which expresses the phase 
non-equilibrium between gas and dust, and p vap the satu- 
ration vapour pressure of graphite. AG is the difference in 
Gibbs free energy between a solid unit and the monomer 
in the gas phase, newly fitted according to the Janaf 
tables (Chase et al. 1985) in units of [J/mol], and R = 
8.31441J mol _1 K _1 the ideal gas constant. The thermal 
b- factors 6* h express departures from thermal equilibrium 
between gas and dust, i. e. Td ^ T g (Gauger et al. 1990, 
Patzer et al. 1998, Woitke 1999). 

Equation @ allows for a generalisation of the defini- 
tion of the thermal stability of dust grains by means of 
t" 1 = 0. This condition states an implicit equation for 
the determination of one of the arguments of r gr , which 
are (utilising Eq.0) Tg,T^,nm\ and ec Consequently, the 
sublimation temperature Tg of a dust grain material can 
be defined as the dust temperature Td which nullifies 
Eq. ©, i.e. r~ 1 (T g , T$, , ec) = 0. In thermal equilib- 
rium (where Td = T g and hence 6* h = 1) , this coincides with 
the usual definition of phase equilibrium (S—l). However, 
in thermal non-equilibrium, the criterion is more compli- 
cate and the resulting sublimation temperatures depend 
on the chemical composition of the gas as well as on the 
considered set of surface reactions. 

2.3. Discrete dust size distribution function 

The dust size distribution function f(a,r,t) is the basis 
for the opacity calculation via Mie theory (see Sect. I2.4|l . 
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Fig. 1. Numerical representation of the discrete size dis- 
tribution function. 

Furthermore, the modelling of the evaporation in the 
frame of the moment method requires the knowledge of 
the dust size distribution function at the lower integration 
boundary f(ag,r,t). We follow the computational method 
developed by Kriiger et al. (1995) to obtain the discrete 
size distribution /, (i= 1, ... , It) on a co-moving grid in 
size space with sampling points a>i(t) (i — 1 , . . . , It + 1 ) , 
where It is the number of grid points at time t, as illus- 
trated in Fig. n Since the growth/shrink rate da/dt is 
independent of a (Dominik et al. 1989), the evolution of 
the size distribution function during the time interval At 
can be obtained by shifting the size grid points uniquely 

by 



Oi(t+At) 

da 
~dt 



, x da 
a lW +- 

«o 



At (i = l,...,/ t +l) 



3r gr (t) 



(13) 



(14) 



The discrete size distribution function is adapted to the 
evolution of nucleation, growth and evaporation according 
to the following numerical scheme 

J(ae,t) > =>• generate new grid point: 

It+At = It + 1 

a It+ i(t+At) = at 
ai t < ag =>■ discard grid point: 

It+At = It — 1 

The value of fj t (see l.h.s. interval in Fig. ^) is obtained 
by integrating the incoming flux of clusters in size space 
over At through the boundary a — ai, assuming that r gr 
and J(ae) are constant within the time step At. Only a 
few more special cases must be considered for this simple 
scheme of a discrete representation of the dust size distri- 
bution function, e. g. when no dust is yet present or when 
the last size bin is to be removed. In these cases, two new 
grid points must be generated/discarded. 

The dust moments can alternatively be calculated by 



KM = — - n -Y % \(ai) i+1 - 



(15) 



which allows for a cross-check with the results obtained 
from the numerical integration of the moment equations, 



where ^ 

f(at,t) = < n n . (16) 
[ , otherwise 

2.4. Opacity Calculations 

The gas is assumed to be optically thin in the con- 
sidered volume and, consequently, gas opacities are ne- 
glected in this model. The dust opacities - in particu- 
lar the dust extinction and angle-integrated scattering co- 
efficients, K^ xt (r,t) and K^ ca (r,i), the albedo 7a (r,t) = 
K A Ca / K A Xt an d the phase function <7a($, r, t) (definition see 
Niccolini et al. 2003) - can be computed by applying Mie 
theory on the basis of the calculated dust size distribu- 
tion function /(a, r, i), if the optical constants of the dust 
material is known (Bohren & Huffman 1983). 

It soon became clear, however, that using exact Mie 
theory in each spatial cell of the model volume (see 
App. A) at every time step according to the calculated 
dust size distribution function would cost too much com- 
puting time. We are therefore forced to use an approx- 
imate procedure instead. Noting that the dust particles 
remain quite small in our model ((a) < 0.1 ^m) the 
small particle limit of Mie theory (Rayleigh limit) can 
be considered, where K^ xt cx K 3 (Gail & Sedlmayr 1987). 
Therefore, as a simplifying approach, we use scaled extinc- 
tion coefficients according to 

*r(r,*) = • (17) 

K% xt and K-j, refer to a standard size distribution, chosen 
to be /(a) oc a -3 ' 5 between 0.001 /im and 1 /im. The advan- 
tage of this procedure is that the Mie calculations must 
be carried out only once per complete model. We use the 
Mie algorithm of Wiscombe (1980) with the optical data 
of amorphous carbon (real and imaginary part of the re- 
fractory index) from Rouleau & Martin (1991) in order 
to calculate K^ xt , 7a and ~g\{&). As a further simplifying 
assumption, the albedo and the phase function are not 
scaled at all, i.e. j\(r,t) = 7 A and g\(<&,r,t) =g\(tf) 5 . 

2.5. Radiative transfer 

The axisymmetric continuum radiative transfer prob- 
lem is solved at each time step by applying our Monte 
Carlo method (Niccolini et al. 2003), which assumes time- 
independence, LTE, angle-dependent coherent scattering 
and radiative equilibrium. 

Several basic improvements of the standard Monte 
Carlo procedure have been introduced in order to reduce 
the Monte Carlo noise in the temperature determination 

5 These relative quantities have little influence on the result- 
ing temperature structure. We have checked that the temper- 
ature differences between taking the true Mie opacities and 
our approximate opacities are less than 5% anywhere in the 
model volume. However, It must be noted that strictly speak- 
ing «;| ca oc Ke in the Rayleigh limit. 



and so to push the computational performance of the 
method. Stellar and cell photon packages are systemat- 
ically generated and independently propagated through 
the dust envelope. The dust particles are assumed to be 
in radiative equilibrium, 

J K f s J A dX = j « A abs 5 A (T d ) d\ , (18) 

where K A bs = «| xt (l — 7a) is the dust absorption coefficient 
and B\ the Planck function. For each cell, the dust pho- 
tons are generated with an initial guess for T^(r,t). Taking 
advantage of the explicit temperature-independence of 
the dust opacities, the correct temperature stratification, 
which satisfies Eq. (|18|l in every cell, is found by iteration 
after the Monte Carlo experiment has been completed, by 
applying a A-operator technique. This high-precision tem- 
perature determination scheme is based on the calculation 
of mean intensities J\ according to the computed path 
lengths of the photon packages through the cells (Lucy 
1999), which gives accurate results for T^(r,t) also in op- 
tically thin situations. For further details about this radia- 
tive transfer method, please consult (Niccolini et al. 2003, 
"method 1"). 

As inner boundary condition for the radiative trans- 
fer problem, a black body sphere with constant radius i?* 
and constant effective temperature T e g is assumed, which 
fixes the bolomctric luminosity — AttR^uT^. A vary- 
ing stellar temperature T+(f) is introduced to determine 
the outgoing flux and the wavelength distribution of the 
stellar photons. As the dust forms and the optical depths 
in the model volume increase, photons that are thermally 
re-emitted or scattered back from the dust shell have a 
certain probability to hit the star, which causes an inward 
directed photon flux at r — R^, (a radiative heating of the 
star) which must be compensated for by an increase of 
the outgoing flux, i.e. T*{t) > T e g. The stellar tempera- 
ture T*(£) is initially set equal to T c g, but is part of the 
aforementioned iteration procedure. Incident radiation at 
the outer boundary is neglected. 

The method is suitable for discontinuous opacity struc- 
tures, such as clumpy media, and is applicable in a broad 
range of optical depths. For the model under consid- 
eration, we use altogether 5 x 10 7 stellar and 5 x 10 7 
cell photon packages on 30 wavelengths points, covering 
0.1 [im ... 250 /im, which results in a mean Monte Carlo 
temperature noise of Amc^cI ~ 0.5 K 6 . One radiative 
transfer calculation takes about 3.5 min on a Cray T3E- 
1200 parallel super-computer using 200 processors. In 
comparison, the integration of the dust moment equations 

6 The Monte Carlo temperature noise Amc^cI has been esti- 
mated by (i) picking a representative ^(restructure from the 
model and smearing it out over angle 9 (see App. A), which 
results in a spherical symmetric dust distribution, (ii) perform- 
ing a full Monte Carlo radiative transfer calculation, (iii) cal- 
culating the standard deviation of the calculated temperatures 
Td(r,6) at constant radii r and (iv) taking the mean value of 
these standard deviations over all radial grid points. 



(Eq. |3J and the calculation ot the size distribution lunc- 
tion (Eq. 1 1 3|> is computationally cheap, consuming only 
about 10% of the total computational time. 

Concerning the determination of the gas temperatures 
T g (r, t), we have no access on detailed molecular opacities 
in our program. Therefore, we simplifyingly assume that 
the gas temperature is given by the black body tempera- 
ture T g {r,t)=T hh (r,t) defined by 

Jj x d\ = jB x (T hh )d\ . (19) 

2.6. Iteration, initial values and time step control 

The model is calculated forward in time by a simple ex- 
plicit iteration scheme. As initial values at t — 0, we choose 
the dust-free case 

f(a,r,t = 0) = & Kj(r,t = 0)=Q. (20) 

The iteration starts by calculating the opacities according 
to the actual values of K 3 (r,t) (see Sect. 12. 4j) . Based on 
these opacities, e.g. K^ Kt (r,t), a complete Monte Carlo 
radiative transfer including temperature iteration is car- 
ried out (see Sect. I2.5|) . which provides the temperature 
structure of the medium at time t: T g (r,t) and Td(r,i). 

Next, the dust moment equations are integrated for- 
ward for a time step At as described in Sect. 12.21 which 
results in the dust moments Kj(r,t + At). This integra- 
tion is performed by means of the assumption that the 
temperatures T g (r,£') and T^(r, t') are constant during 
t' e [t,t + At]. 

The moment equations are numerically integrated sep- 
arately for each cell (see App. A). An internal time step 
control is necessary to account for rapid changes of the 
dust quantities during At and the stiff behaviour of the 
dust component close to the sublimation temperature. 
The carbon abundance ec(r,t'), the concentrations of the 
chemical species rifc(r, t') and the discrete size distribution 
f(a,r,t') are essential parts of this integration. Once the 
dust-chemistry time-step is completed, the model proceeds 
with a new opacity calculation (see above), providing the 
basis for the next radiative transfer and so on. The model 
is calculated forward until the spatial dust structures have 
fully developed. Typically, this occurs after 150 — 300 time- 
steps, which requires altogether 2 — 3 eight-hour-jobs on 
the Cray T3E-1200 parallel super-computer using 200 pro- 
cessors. 

The choice of the time step At is adapted to the phys- 
ical situation at time t in order to assure that T g (r,t') 
and Td(r, t') in fact remain approximately constant dur- 
ing t' G [t,t+At\. We achieve this goal by combining sev- 
eral time-scale criteria. Most importantly, K3 must not 
change too much during At in any cell, which would re- 
sult in a relevant opacity-change that might influence the 
temperature structure. The maximum allowed change of 
K3 during At is limited by introducing an absolute toler- 
ance (2-10~ 3 [e^ — eo]) an d a relative tolerance (5%). 



j. Kesults 

The results obtained from our axisymmetric model cal- 
culations are presented in the following way. We start by 
showing one exemplary model and explain the main phys- 
ical processes in Sect. 13. ll The radial development of the 
dust is further discussed in Sect. 13.21 and the angular devi- 
ations from it (structure formation) in Sect. 13.31 Sect. 13.41 
gives more insight into the nature of the radiative/thermal 
instability and discusses where the dust forming gas in an 
AGB star wind might be affected by this instability. The 
functional dependencies of the introduced parameters are 
outlined in Sect. l3.5l and Sect . 13 . fil discusses the influence of 
the forming dust clumps on the spectral appearance of the 
star by means of calculated spectral energy distributions 
and monochromatic images. 

3.1. One exemplary model 

The development of the forming dust shell in a typical 
model, in terms of the degree of condensation f C ond(r,t) 
and the dust temperature Td(r,t), is shown in Fig.|3 
The denser regions close to the star condense first, be- 
cause the nucleation and growth of the dust particles pro- 
ceeds faster at larger densities. Consequently, the spatial 
dust distribution, at first, resembles the slightly inhomo- 
geneous gas distribution in the circumstellar shell with a 
cutoff at the inner edge as a consequence of the tempera- 
tures being too high for nucleation close to the star, and 
a smooth outer boundary due to the radially decreasing 
density. According to the different choice of the density 
inhomogeneities above/below mid-plane (see Ea. lA.6(l . the 
resulting spatial variations in the degree of condensation, 
Ag log /cond, are larger above than below the mid-plane at 
early phases of the model (see upper left plot in Fig.[2j). 
Noteworthy, Ag\ogf con d is larger than Aglogn^), be- 
cause the inverse dust formation time-scale initially scales 
with the density to a power of 3 ... 4 (Woitke 2001). 

The process of dust formation continues for a while 
in this way, until the first cells close to the star become 
optically thick. Each already optically thick cell casts a 
shadow into the circumstellar envelope wherein the tem- 
peratures decrease by several 100 K which improves the 
conditions for subsequent dust formation therein. At the 
same time, scattering and re-emission from the already 
condensed cells intensifies the radiation field inbetween. 
The stellar flux finally escapes preferentially through the 
segments which are still optically thin, thereby heating 
them up and worsening the conditions for further dust 
formation there. These two contrary effects amplify the 
initial spatial contrast of the degree of condensation intro- 
duced by the assumed density inhomogeneities. In the end, 
radially aligned, cool, linear dust structures, henceforth 
called dust fingers have developed, which point towards 
the star and are surrounded by warmer, almost dust-free 
regions at the inner edge of the forming dust shell. The 
lengths of the fingers is of the order of 0.5 R±. 
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Fig. 2. Self-organisation and spatio-temporal structure formation, triggered by radiative/ thermal instabilities, during 
the formation of a dust shell around a carbon-rich AGB star. The figures show contour plots of the degree of conden- 
sation log(/cond) (left column) and the dust temperature structure Td (right column) in the x/z-plane (cut through 
the y = 0-plane) at three selected time steps as indicated on top. On the r.h.s., black colour indicates regions that are 
not included in the model. The white circles in the centre of each figure mark the star. Because of axisymmetry, all 
data points appear twice at ±x. Parameters: T e g = 3600K, L^ — 3000 Lq, i?* = 9.79 • 10 12 cm, C/O — 2, H p = 0.25 i?*, 
™<H}M = 3.7-10 10 cm- 3 and C* 2 = 0.1/0.03. 



6.Z. I he chemical wave 

Apart from the self-organisation of the dust in the angular 
direction, the model reveals basically the formation and 
evolution of a radial dust shell. For a better visualisation 
of these processes, we have plotted several angle-averaged 
quantities (X)g(r) in Fig. The effective formation of 
dust generally requires a suitable combination of gas den- 
sity and temperature, called the dust formation window 
(Gail & Sedlmayr 1998). In the beginning of the model, 
favourable temperature conditions for efficient nucleation 
are only present in a restricted radial zone close to the 
star. However, as time passes, the dust shell becomes opti- 
cally thick which dams the outflowing radiation and leads 
to a considerable increase of the temperatures inside of 
the shell (radiative backwarming) . Consequently, the zone 
of effective dust formation shifts outward with increasing 
time. Moreover, the temperatures at the inner edge of the 
shell temporarily exceed the sublimation temperatures T$ , 
and the dust shell begins to re-evaporate from the inside. 
The upper plot of Fig. [3] shows that Td in fact temporar- 
ily exceeds its equilibrium value which is reached at later 
stages (see Sect. 13. 4J) . 

These two effects result in an apparent motion of the 
dust shell, driven by dust formation at the outer edge 
and dust evaporation at the inner edge of the shell, 
which will be denoted by chemical wave in the follow- 
ing 7 . Chemical waves are a well-known phenomenon in the 
laboratory, for example the famous Belousov-Zhabotinsky 
reaction (Zaikin & Zhabotinsky 1970), when chemical sys- 
tems show oscillatory solutions, sometimes radiatively 
controlled (Schebesch & Engel 1999), even without fluid 
bulk velocity fields. The apparent propagation velocity of 
the chemical wave, calculated from Fig. decreases from 
w 2km/s (early phase) to w O.lkm/s (late phase). 

Figure 0] shows more details of the structure of the 
chemical wave. The outer regions ahead of the wave are 
featured by a positive nucleation rate J(ai) = J* > and 
a positive growth time-scale T gr > 0. In the centre of the 
wave, /cond is maximum and J(ai) and 1/T gr vanish. The 
layers behind the wave (close to the star) usually pos- 
sess a negative flux through the size integration boundary 
J(ai)<0 and a negative growth time-scale r gr < 0. 

After the passage of the chemical wave, the system re- 
laxes towards a stable equilibrium state, where the dust 
and the temperature structure of the medium are mutu- 
ally coupled in a fine-tuned way. The equilibrium is en- 
forced by a self -regulation mechanism: More dust causes 
an increase of the temperatures in the neighbourhood via 
backwarming, which leads again to dust evaporation and 
vice versa. This mechanism results in a physical state close 
to phase equilibrium T c i — > Ts (or, alternatively speaking, 
r gr — > oo, see Sect. 12.2(1 for long times after the chemi- 
cal wave has passed, featured by a considerable increase 

7 In order to avoid confusion, we want to stress once more 
that no bulk velocities exist in our model (v = 0). We will 
nevertheless use the technical term "chemical wave" to denote 
wave-like propagations of chemical fronts. 
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Fig. 3. The chemical wave: formation and propagation of 
a dust shell via dust formation at the outer edge and 
dust evaporation at the inner edge. The figure shows 
the time evolution of the angle-averaged dust temper- 
ature (Td)e(r) (upper plot) and degree of condensation 
(/cond)e(^) (lower plot). The quantities are shown for it- 
eration step it = (t = 0, full), it = 10 (t = 0.21 yr, dot- 
ted), it = 20 (t = 0.36 yr, dashed), it = 40 (t = 0.84 yr, 
dashed-dotted) , it = 70 (t = 2.1 yr, dashed-triple-dotted), 
it = 120 (t = 5.7 yr, long-dashed) and it = 195 (t = 13 yr, 
full). 



of (f C ond)e( r ) an d a moderate decrease of {T g )g(r) with 
increasing r. The relaxation time-scale towards this equi- 
librium state is approximately given by the characteris- 
tic dust growth/evaporation time-scale, \K 3 / (dK 3 /dt)\, 
which is density-dependent. Limited by our explicit nu- 
merical iteration scheme, we are forced to keep the compu- 
tational time-step At smaller than this characteristic time- 
scale in the inner regions (see Sect. I2.6|) which makes it 
computationally difficult to trace the chemical wave much 
further out, where the dust formation proceeds slower by 
orders of magnitude due to the lower densities. 
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Fig. 4. Detailed radial snapshots of various angle-averaged quantities. Upper boxes: dust temperature (Td)e (full 
line) and gas temperature (T g )g (dashed). Second boxes: positive nucleation rate (J*/«(H))e (full line) and negative 
evaporation rate —(J(at)/n(n\}g (dashed). Third boxes: inverse growth time scale (l/r gr )g. Lower boxes: degree of 
condensation (/ CO nd)e- The errorbars indicate no errors but standard deviations AqX of the angular variations of the 
calculated quantities X(r, t). 



3.3. Structure formation 

Beside the evolution of the dust component in the 
r-direction, represented by angle-averaged quantities 
(X)g(r), Fig. 01 shows the standard variations of angu- 
lar variations of the calculated quantities by means of er- 
rorbars. These "errorbars" have nothing to do with real 
"errors" 8 , but give an impression of the magnitude of the 
angular variations occurring in the model (compare also 
Fig. 0). 

In general, the chemical wave is found to leave behind 
a strongly inhomogeneous dust distribution, where / con d 
varies between zero and particular maximum values, which 
depend on the gas density and the distance to the star. 
From Fig.|H(r.h.s.) we see that usually A$X(r) > (X)$(r), 
where X G { J(ai), 1/V gr , /cond}) f° r long times after the 

8 In fact, statistical errors in the temperature determination, 
AMcTd Ss 0.5 K, occur in the model due to the Monte Carlo 
noise (see Sect. 12.51 which are, however, much smaller than 
the depicted angular variations AgT. 



passage of the chemical wave in the regions close to the 
star, i. e. the angular variations of all relevant dust quanti- 
ties are usually larger than their mean values. This result 
is a complicate consequence of several factors. 

1. The dust formation process requires the formation of 
seed particles (nucleation), which depends exponen- 
tially on T g with threshold character. Therefore, small 
temperature variations may result in large variations 
of the dust quantities. 

2. The dust growth rate is density-dependent. Small den- 
sity variations temporarily lead to large contrasts of 
/cond in early phases of the condensation process. 

3. The sublimation temperature Ts depends on the den- 
sity. Denser regions are more resistant against thermal 
evaporation, because the partial pressures of carbon 
molecules in the gas phase are higher and hence the 
supersaturation ratio S larger for larger gas densities. 



4. Kadiative transfer effects introduce a non-local spa- 
tial coupling. In particular, dense cells with high / con d 
cast shadows. In the shadowed regions behind these 
"clouds", the conditions for the subsequent formation 
(nucleation) and the survival of the dust (thermal sta- 
bility) are improved due to shielding effects. 

The first two points lead to an inhomogcneous dust dis- 
tribution already during early phases of the dust con- 
densation process. As soon as the dust shell becomes 
optically thick, the backwarming enforces a partial re- 
evaporation of the dust, until the radiative energy finds 
ways to partly escape the system. Due to the threshold- 
like temperature-dependence of the dust survival, and the 
density-dependence of Tg, the dust in the slightly denser 
cells close to the star just even survives, whereas the dust 
in the adjacent (slightly thinner) cells evaporates com- 
pletely. Figure 0] shows that 1 /r gr (zero indicates phase 
equilibrium) varies from slightly positive to considerable 
negative values behind the chemical wave, i. e. these re- 
gions consist of a mixture of cells where the dust is just 
even thermally stable and cells where all dust particles 
evaporate/have evaporated. 

The spatial coupling via radiative transfer effects 
brings order into this chaos. If the dust in one particular 
cell can resist the radiation field, it shields the radiation 
and facilitates the survival of the dust in the shadow of this 
cell. Therefore, the final equilibrium state of the medium 
can in fact possess a non-trivial spatial structure, where 
cool, optically thick segments can coexist besides warmer, 
optically thin segments. The radiative flux finally escapes 
preferentially through these optically thin segments where 
the degree of condensation remains low (see upper right 
plot of Fig. 0J. Thus, finger- like dust structures develop 
(see Fig. |2J which point toward the star. 

3.4. The metastable region 

The driving mechanism for the structure formation is an 
unstable physical control loop as thoroughly discussed in 
Paper I and once again sketched in the upper part of 
Fig-El The physical feedbacks [1] and [3] (temperature- 
dependence on radiation field and opacity-increase in case 
of condensation) are always positive and active. However, 
the feedback [2] diminishes at too low temperatures: Let 

° o 

/condM and T(r) denote the degree of condensation and 
the dust temperature structure in the final equilibrium 
state reached long after the passage of the chemical wave. 
The ability of the medium to condense can then be ex- 

o 

pressed by (f C ond)(T). At too low temperatures (large r), 
a small temperature change 5T, e. g. by shadow casting, 
has only little influence on this quantity (as indicated by 
the vertical arrows in Fig.EJ), the control loop is broken 
and an outer boundary of the metastable region results 
(see lower part of Fig.[5l). On the other hand side, the 
feedback [4] (the non-local effect of opacity on the tem- 
perature structure, in particular shadow casting) requires 
that the medium is not completely optically thin, which 
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Fig. 5. Sketch of the radiative/ thermal instability of dust 
formation and related size of the metastable region. 

creates an inner boundary of the metastable region. Here, 
the achievable degree of condensation is too small because 
of too high temperatures, and since gas opacities are unim- 
portant, the medium cannot attain a considerable optical 
thickness. 

Consequently, significant structure formation only oc- 
curs in a limited spatial zone of the model (see Fig. [2j 
between two boundaries, denoted by [4] and [2] in Fig. [5] 
henceforth called the metastable region. Here, the matter 
finally remains in a metastable state which is not fully con- 
densed nor completely dust-free, where small fluctuations 
cause large effects. 

From the numerical results, we infer that the size of 
this metastable region is essential for the significance of 
the forming dust structures. Only if this region is large, 
there is enough space for the radiative/thermal instabil- 
ity to bear in fact large-scale spatio-temporal structures. 
This size is related to the radial density gradient of the 
medium. With outward decreasing temperature, a certain 
radial gradient of the vapour pressure of the dust grain 
material Pvap (2d (?")) is given 9 . If the gas density decreases 
with a similar gradient, the supersaturation ratio S is of 
the order of unity in an extended radial zone. 

3.5. Functional dependence on model parameters 

This paragraph describes the influence of the model pa- 
rameters on the structure formation process, based on the 



9 Noteworthy, the Td-gradient is much smaller inside of an 
optically thick dust shell than in the optically thin case. 



experience gained irom other model calculations with dif- 
ferent parameters not shown in detail in this paper 10 . The 
numbers in squared brackets in the item list below are the 
parameter values of the standard model discussed so far 
(in Figs. 02 and EJ| . 

T e g [3600 K]: The effective temperature of the star trig- 
gers the mean radial distance of the dust formation zone. 
If T c ff is decreased, the metastable zone (see Sect. 13 .4|) is 
located closer to the star and results to be narrower, which 
leads to the formation of less significant dust structures. 

fi'(H) (fo) [3.7 • 10 10 cm" 3 ]: The gas density in the circum- 
stellar environment is responsible for the dust growth 
time-scale and, hence, for the velocity of the chemical 
wave. By changing this parameter, the model mainly pro- 
ceeds on a different time-scale. Larger densities also lead 
to an increase of Ts, which means that the formation 
and the survival of dust is already possible closer to the 
star. Furthermore, larger densities potentially lead to the 
formation of more dust on an absolute scale, which pro- 
duces larger optical depths. In models with strongly re- 
duced n(H>(' , o) ) the resulting optical depths are too small 
to cause significant temperature changes via backwarm- 
ing and shadow formation, and the system decouples 
spatially, i. e. every cell condenses on its own, without 
much feedback on the neighbouring cells. In models with 
larger ri( H ^(ro), we observe that only slightly larger opti- 
cal depths occur as compared to the standard model, but 
that the dust shell is narrower. This is a consequence of 
the re-evaporation process which sets in as soon as a cer- 
tain critical value of Ti Mm is reached (here about 1...3). 
The temperatures temporarily reach higher values behind 
the wave as depicted in Fig. and the dust evaporation 
is faster and more complete on the inside of the shell. In 
both cases of varied n(H}( r o), less pronounced structure 
formation is found to occur. 

C/O [2.0]: The initial carbon-to-oxygen ratio is a mea- 
sure for the amount of condensable material in the gas. 
Consequently, changing C/O has a similar influence than 
changing the overall density level by means of n^(ro). 

Hp [0.25 R*]: The scale height is important for the ra- 
dial extension of the metastable zone (see Sect. 13.4]) . We 
observe from models with larger H p that the size of the 
metastable zone is smaller and, consequently, the struc- 
ture formation results to be less significant. However, even 
for models with ocl/r 2 , where the density gradient is 
much shallower, some structures still develop. 

C2 [0.1/0.03]: The degree of density inhomogeneities (see 
Eq. IA.4fl is not a critical parameter of the model ! From 
Fig. [21 we see that the number and the shape of the dust 
structures is finally similar above and below the mid- 
plane, despite the different values for C2. The finger-like 
structures even develop for C2 = 0. For larger density in- 

10 See mpeg-movies at http://astro.physik.tu-berlin.de/ 
~woitke/sfb555 . html 



homogeneities, the structures become slightly longer and 
less numerous. 

Spatial resolution [71 x 71 grid points]: The number of ra- 
dial grid points in the model is important for the proper 
resolution of the radiative transfer and the propagation 
velocity of the chemical wave. The number of angular grid 
points is crucial with respect to the shape and, in par- 
ticular, the width of the resulting dust structures. In the 
depicted model (Fig. EJ), -ZV struc = 27 dust structures can 
be identified within Ng + 1 = 71 angular grid points, i. e. 
the resulting number and width of the dust fingers is at 
the limit of the angular resolution of the model. 

The resulting width of the dust fingers is possibly a conse- 
quence of the way how we have introduced the density in- 
homogeneities in our model. Using random numbers for all 
cells (see App. A), the resulting density inhomogeneities 
are spatially uncorrelated and the characteristic scale of 
the density variations is naturally given by the size of the 
cells, which is apparently important for the width of the 
dust fingers. We have run model calculations with up to 
181 angular grid points (which forces us to reduce N r to 
40 due to computer memory and time-consumption limi- 
tations), where the finger- like dust structures become nar- 
rower at the head and slightly conic according to the shape 
of the penumbra. The ratio N stvU c/Ng becomes less if N$ is 
increased. However, definite predictions about the diame- 
ter of the dust structures are difficult to make, because of 
the limited angular resolution in our model. 

3.6. Spectral appearance 

From an observational point of view, it is interesting 
to discuss how a "clumpy" dust distribution around an 
AGB star would appear both spectroscopically and in 
monochromatic images. Figure shows the calculated 
spectral energy distribution (SED) of the standard model 
at the end of the computation. The depicted flux has 
been averaged over all escape directions 9 CSCl 4> csc (see 
Niccolini et al. 2003). The resulting SED is featureless ac- 
cording to the amorphous carbon opacities and can rea- 
sonably well be fitted by a black body with colour tem- 
perature k 2700 K ... 2900 K (note the difference to the 
effective temperature T e g = 3600K). 

The Monte Carlo method allows for an identification 
of the kind of photons which reach the observer. The op- 
tical (blue) part of the spectrum mainly results from di- 
rect star light, which is attenuated by dust extinction, and 
from multiply scattered star light. The near IR spectral 
region around 2 /im ... 4 /im is dominated by thermal dust 
emission, whereas for longer wavelengths, the dust enve- 
lope becomes optically thin and the direct star light again 
becomes the most important contributor 11 . 

To our surprise, the resulting angle-averaged SED con- 
tains almost no evidence for the clumpy dust distribution 

11 This is an artifact of our model, because it is radially not 
sufficiently extended. 
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Fig. 6. Calculated spectral energy distribution of the 
clumpy model at t = 13.2 yr, averaged over all escape direc- 
tions. Flux values correspond to a distance of 1 R+. Full 
line: calculated total flux, dotted line: direct star light, 
thin dashed line: scattered star light, dashed-dotted 
line: thermal dust emission (partly scattered), thick grey 
dashed line: total flux of a spherically symmetric refer- 
ence model (see text). 

Table 2. Comparison of the standard "clumpy" model to 
a spherically symmetric reference model at i = 13.2yr (see 
text). (Td)M d is the mass-averaged dust temperature and 
Tipj, the radial optical depths at A = l/im. 





clumpy model 


reference model 


M gas [M Q ] 


1.35 • 10 -6 


1.35 • 10 -6 


Af dust [Mq] 


5.23 • 10" 10 


4.69 ■ 10 -10 




0.51 ... 2.7 


0.98 


T*[K] 


3610.7 


3610.6 


<Td>M d [K] 


1638 


1654 



in the model, contrary to what results from simplified 
models for spectral analysis (e. g. J0rgensen et al. 2000) . 
We have calculated a time-dependent reference model 
which is identical to the standard model, except for a 
choice of Ng = 1, i.e. all "cells" are closed spherical shells 
in the reference model, which enforces the results to be 
spherically symmetric. Figure shows that the calculated 
SED of the "clumpy" standard model is practically indis- 
tinguishable from the SED of this reference model at an 
comparable instant of time. There is only a slight excess in 
the blue part of the spectrum as compared to the reference 
model (e.g. +15% at A = 0.34/im, +4% at A = 0.55 /mi, 
+ 1% at A = 0.72 /mi), caused by the increased average es- 
cape probability in the blue in case of a clumpy medium. 

Tabic [21 shows more details about this comparison. In 
fact, about 10% more dust (by mass) condenses in the 
standard model as compared to the spherically symmet- 
ric reference model. The standard deviation of the op- 
tical depth at A = 1 /im is larger than its mean value 



\Tif_ini). Ine mass-averaged dust temperature {ld/M d — 

(E e ^nf H) ^|T d «)/(E e n^H)^l) ^ slightly lower 
than in the reference model, where is the volume of 
cell £. In summary, if the assumption of spherical sym- 
metry is relaxed to axisymmetry, more dust condenses in 
the stellar environment which forms and survives in ra- 
diatively shielded, cool clumps, resulting in a lower mean 
dust temperature. However, the characteristic changes are 
only of the order of 10%, and it is noteworthy that the dif- 
ferent single effects tend to compensate each other in view 
of the spectral appearance. 

From Fig. we can infer that the spectral region 
around 2 /im ... 4 /im, which approximately coincides with 
the maximum of the thermal emission from the dust en- 
velope, is most promising to reveal the spatial dust distri- 
bution around the star by direct imaging. Figure shows 
simulated pictures at A = 2.2 /mi. The light from the dust 
shell mainly arises from the glowing heads of the dust fin- 
gers, which appear as rings around the star due to our 
assumption of axisymmetry. The different physical contri- 
butions are shown separately in Fig. [7] However, the in- 
tensities received from the dust shell are about one order 
of magnitude less than the intensities received from the 
star. The lower right image depicts the total brightness 
distribution, where basically only a stellar disk is visible, 
which is partly obscured by the inhomogeneous dust dis- 
tribution around the star. 

Observations of such non-uniformly bright stellar disks 
are usually explained by cool and/or hot spots on the 
stellar surface, e. g. in case of the cool supergiant a Ori 
(Lynds ct al. 1976, Gilliland & Dupree 1996). However, an 
inhomogeneous dust distribution around the star (which 
is a perfectly uniform black body sphere in our model) 
can obviously lead to a similar optical appearance, even 
if there is not much direct evidence for dust in the im- 
age 12 . Thus, the interpretation of monochromatic images, 
e. g. by speckle interferometry, via cool and hot spots on 
the stellar surface is delicate if there is only the slightest 
evidence for the existence of circumstellar matter, e. g. by 
spectroscopy, which could be inhomogeneous. 

In fact, the direct + scattered star light only provides 
« 40% + 8% of the total signal integrated over the en- 
tire image, whereas ~ 37% + 15% originates from the 
dust envelope (direct + scattered dust emission). However, 
the dust emission is spread over a much more extended 
area, here ~ 7r(4i?*) 2 , than the direct star light ~ 7ri? 2 . 
Consequently, the intensities received from the dust en- 
velope are generally much weaker than the intensities re- 
ceived from the stellar disk (factor ~16). 



On the basis of only one observational image, it seems ac- 
tually difficult to decide which physical interpretation is cor- 
rect (surface convection/non-uniform dust formation). Multi- 
epoch, or even better simultaneous multi-wavelength images 
would be required to come to definite conclusions, since sur- 
face convection and dust formation may occur on similar time- 
scales. 
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Fig. 7. Simulated images at 2.2 /jm for an inclination angle of 20° between the equator and the observer's direction 
of the innermost ±4.5 i?* of the model. The dust distribution refers to the final state of the model [t = 18.4 yr). The 
total signal is shown in the lower right plot, whereas the other (artificial) images show the different constituents of 
the signal (total = direct star light + scattered star light + direct dust emission + scattered dust emission). Note the 
different scaling of the intensities, given in units of 10 Wm _2 /im _1 for a reference distance of 1 i?*. The picture of the 
direct star light is very similar to the total signal depicted on the lower rights, and is hence omitted. 



Ihe different physical contributions to the image are 
separately shown in Fig. [7| A composite picture of scat- 
tered star light + direct dust emission + scattered dust 
emission (lower left + upper left + upper right in FigQ) 
would simulate a coronographic image where the stel- 
lar disk is blinded out. Such observations could in fact 
reveal the spatial dust distribution around the star. Of 
course, since our model is axisymmetric, the "clumpy" 
dust distribution in the model appears as superposition of 
rings in the image, which is not expected to be realistic. 
Nevertheless, the images give a first impression of the kind 
of brightness distributions to be expected from a cloudy 
circumstellar environment. 

4. Conclusions and discussion 

This paper has shown that the formation of dust shells 
in the circumstellar environments of late-type stars can 
be unstable. Spontaneous symmetry breaking may occur 
due to a radiative/thermal instability in the dust form- 
ing gas, which leads to the development of cloud- like dust 
structures close to the star. 

These results have been obtained on the basis of time- 
dependent axisymmetric models, which combine a kinetic 
description of carbon dust formation/evaporation with de- 
tailed, frequency-dependent radiative transfer by means of 
a Monte Carlo method, in the static case. The simulations 
show that the dust preferentially forms behind already 
condensed regions, which shield the stellar radiation. In 
the shadow of these clumps, the temperatures are lower 
by a few 100 K which triggers the subsequent formation 
and facilitates the survival of the dust close to the star. As 
final result, numerous finger-like dust structures develop 
which may have an radial extension as large as 0.5 i?* and 
point towards the centre of the radiant emission, similar to 
the "comctary knots" observed in planetary nebulae and 
star formation regions. 

The cloudy dust distribution has little effect on the 
calculated spectral energy distribution of the star, in com- 
parison to a spherically symmetric model, but signifi- 
cantly influences the optical appearance of the circumstel- 
lar environment in near IR monochromatic images (e. g. at 
A = 2.2/im). In particular, an inhomogeneous dust distri- 
bution around the star leads to a likewise non-uniformly 
bright appearance of the stellar disk. Comparable obser- 
vations of late-type giants are commonly interpreted in 
terms of hot/cool spots on the stellar surface. However, 
our model suggests that a different physical explanation 
by dust is possible for these observations, even if the dust 
is barely visible in the image. 

Due to computational time constraints, we have so far 
only been able to study the static case where velocity 
fields are ignored. In this case, the main feature of the 
model is the formation of a chemical wave which radi- 
ally propagates outward, driven by dust formation on the 
outer edge and dust evaporation due to backwarming at 
the inner edge. The optical depth of this chemical wave 
reaches about Ti^ m ^ 1 ... 3 with a density-dependent prop- 



agation velocity ol « 0.1 km/s ... I km/s. Ihe wave leaves 
behind a strongly inhomogeneous dust distribution close 
to the star, where the medium relaxes towards phase equi- 
librium. However, this process is unstable and results in 
the aforementioned simultaneous occurrence of cool dusty 
(optically thick) segments beside warmer, almost dust-free 
(optically thin) segments, through which the radiative flux 
finally escapes preferentially. 

According to dynamical (but spherically symmetric) 
models of dust-forming AGB stars (e. g. Winters et al. 
2000, Schirrmacher et al. 2003, Homer et al. 2003, Sandin 
& Hofner 2003), the formation of dust mainly occurs in 
particular phases of the model triggered by the pulsation 
of the star, which results in radial dust shells. During such 
a dust shell formation phase, the effect of a re-evaporation 
from the inside is a typical feature, similar to the be- 
haviour of our chemical wave. According to the present 
paper, this process should be unstable and might result in 
departures from spherical symmetry. 

We believe that this instability can provide a basis for a 
better understanding of inhomogeneous dust distributions 
showing up in many observations. However, direct predic- 
tions for particular objects are difficult to make, because 
hydrodynamics is missing in the current model. Since the 
dust is blown away as soon as it forms, the system has only 
little time to relax towards phase equilibrium at the inner 
edge of the dust shell. On the one hand, this time may be 
too short to produce well-grown spatial dust structures as 
shown in this paper. On the other hand, even small de- 
viations from spherical symmetry may trigger important 
dynamical effects, e. g. 

— The radiation pressure on dust grains, which is deliv- 
ered to the gas via frictional forces, mainly depends 
on the radiative flux and the degree of condensation. 
If slightly more condensed regions exist, they might be 
accelerated outward, whereas less condensed regions 
stay behind or even fall back. 

— Optically thick dust clouds will be confined by radi- 
ation pressure, because the bolometric radiative flux 
is larger at the inner edge of the cloud facing the star 
than at its self-shielded outer edge (this effect does not 
occur in spherically symmetric models where r 2 F(r) = 
const in radiative equilibrium). Driven, pancake-like 
structures could evolve due to this effect. 

— The hydrodynamical process of cloud acceleration due 
to radiation pressure is not well-studied, apart from the 
special case of spherical symmetry. This process may 
be dynamically unstable itself (e. g. Ray leigh- Taylor, 
Kelvin-Helmholtz). Velocity disturbances generated by 
these instabilities may have an important feedback on 
the dust forming medium. 

Thus, we propose a new hypothetical scenario for dust- 
driven AGB star winds: Excited by hydrodynamical, ra- 
diative or thermal instabilities, dust clouds are formed 
from time to time close to the star in temporarily shielded 
areas, which are accelerated outward by radiation pres- 
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Fig. A.l. Spatial grid of the axisymmetric model. Each 
cell is a torus which appears twice (at x > and x < 0) 
in this cut through the y = plane. Blank regions are not 
included in the model volume (the photons are assumed 
to pass through these regions). The star is indicated as 
central disc. 



sure. At the same time, thinner, dust-free matter falls back 
towards the star at different places. A highly dynamical 
and turbulent environment close to the star would be cre- 
ated in this way, which can be expected to bear again a 
strongly inhomogeneous dust distribution. 

In order to verify this hypothetical scenario, much 
more elaborate model calculations would be required 
which may well exceed the present capabilities of parallel 
super-computers. Various processes must be traced in 3D, 
using dynamical models with detailed radiative transfer 
and time-dependent dust-chemistry 
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Appendix A: Spatial grid 

The model volume is subdivided into fixed spatial cells, 
where the cell boundaries are specified in spherical coor- 
dinates r, 9 and </>, using prescribed sampling points for 
the radius and the latitude angle, r n and 0k, respectively. 
A cell is the volume defined by the following coordinate 



r e [r n -i,r n \ 

9 6 [Ok-l,0k] 

4> e [0,2tt] . 



n = 1, 
k = 1, 



(AT) 
(A.2) 
(A.3) 



According to Eq. (|A.3I) the cells are closed tori, i. e. the 
model is axisymmetric (2D). Each cell is specified by one 
radial and one angular index, n and k, respectively, which 
are abbreviated by a multi- index £ = (n, k). For the model 
discussed in detail in this paper, we choose a radial grid 
with N r + 1 = 71 logarithmic equidistant sampling points 
between r /i?* = 3.0 and rjy r /-R* = 5.5 and Ng + 1 = 71 
equidistant angular grid points between 0q= and 0^ g = 
it. The spatial grid is visualised in Fig. IA.1I 

The prescription of the density structure in the model 
volume (Eq. ^) is technically realised in each cell £ via 



(H) 



C\ ■ exp 



C 2 Z 



(A.4) 



where r ' = (r„_i+r„)/2 the mean radial distance of cell £ 
and Z a random number with Gaussian distribution and 
unity variance 



v/-21n(l - z{) ■ sin(27r z 2 ) . 



(A.5) 



z i, z i £ [0, 1) are equally distributed pseudo random. C\ = 
6 • 10 15 cm -3 is a constant chosen to produce the desired 
mean hydrogen nuclei particle density at the innermost 
radial grid point j^h^o) 13 - In order to study the influence 
of the assumed level of density inhomogeneities, we choose 
the second constant as 



C 2 = 



0.10 , above mid-plane 
0.03 , below mid-plane 



(A.6) 



The standard deviation of the density at constant radius, 
A g log n ( H), is 10%doge«4.3% above and 3%doge«1.3% 
below the mid-plane of the model volume (z = 0). 
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